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This article describes a purely analytic approach to urn models 
of the generalized or extended Polya-Eggenberger type, in the case 
of two types of balls and constant "balance," that is, constant row 
sum. The treatment starts from a quasilinear first-order partial dif- 
ferential equation associated with a combinatorial renormalization 
of the model and bases itself on elementary conformal mapping ar- 
guments coupled with singularity analysis techniques. Probabilistic 
consequences in the case of "subtractive" urns are new representa- 
tions for the probability distribution of the urn's composition at any 
time n, structural information on the shape of moments of all orders, 
estimates of the speed of convergence to the Gaussian limit and an ex- 
plicit determination of the associated large deviation function. In the 
general case, analytic solutions involve Abelian integrals over the Fer- 
mat curve x +y h = 1. Several urn models, including a classical one 
associated with balanced trees (2-3 trees and fringe-balanced search 
trees) and related to a previous study of Panholzer and Prodinger, 
as well as all urns of balance 1 or 2 and a sporadic urn of balance 3, 
are shown to admit of explicit representations in terms of Weierstrafi 
elliptic functions: these elliptic models appear precisely to correspond 
to regular tessellations of the Euclidean plane. 

0. Introduction. In this study, we revisit the most basic urn model, 
namely the "generalized" (or "extended") Polya-Eggenberger urn model 
with two types of balls, as described in the reference book of Johnson and Kotz 
(1977). Under this model an urn may contain two types of balls, say "black" 
(B) and "white" (W). The composition of the urn at time is fixed. At time 
n, a ball in the urn is randomly chosen and its color is observed (thus the 
ball is selected, examined and then placed back into the urn): if it is black, 
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then q black and (3 white balls are subsequently inserted; if it is white, 
then 7 black balls and 5 white balls are inserted. The evolution rule is then 
summarized by a 2 x 2-matrix 



drawn 


added 


I 


B W 


B 


a (3 


W 


7 5 



Negative values of the diagonal entries a, 5 are permissible and interpreted 
as an extraction (rather than an insertion) of balls; a model with both di- 
agonal entries negative will be called here an urn with subtraction (of balls 
of the color chosen). The off-diagonal entries /3, 7 are always taken to be 
nonnegative. 

The urn model is said to be balanced if a + (3 = 7 + <5, in which case 
the common sum of the matrix rows is the balance, denoted throughout 
by s. The 2x2 urn model may lead to widely differing behaviors de- 
pending on the values of the integer entries ct,/3, 7 and 5. For instance, 
Kotz, Mahmoud and Robert (2000) mention the (balanced) urn with ma- 
trix (g j) for which the number of white balls picked in n steps grows 
stochastically like n 1//4 . Strikingly, Kotz, Mahmoud and Robert (2000) also 
study the (unbalanced) urn associated to (\ J) and show the correspond- 
ing number to be ~n/logn in probability under a Poisson model. We do 
not address in this article models with more than two colors; see the paper 
of Smythe (1996) for a thorough probabilistic treatment, the works of Aldous 
(1991), Aldous, Flannery and Palacios (1988) for a discussion of almost sure 
convergence issues, and the comprehensive and independent recent studies 
of Janson (2004, 2005). 

Our interest throughout this article is in urn models that are balanced. 
The conditions of having a matrix 

(1) M= (" 5) with a + /? = 7 + 5 = s,/3>0,7>0, 

are invariably assumed. We also allow ourselves on occasion to describe M 
linearly as (a,(3;j,5). In such a case, each elementary action on the urn 
results in having the total number of balls increase by the fixed quantity s, 
so that the population at time n has a predictable cardinality, which is ex- 
actly to + sn if to is the initial size at time 0. For urns involving subtraction, 
certain simple arithmetic conditions on the parameters, called tenability (the 
Webster dictionary defines "tenable" as meaning "capable of being main- 
tained" ) , ensure that the process cannot be "blocked" ; these conditions are 
recalled in Section 1, (6), and are assumed to hold. 

Balanced 2x2 urn models have been in particular considered by Bagchi and Pal 
(1985) who show the following: under a supplementary technical condition, 
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namely that the ratio between eigenvalues of the matrix, (a — (3)/(a + f3), 
lies in (—00, k), the distribution of the number of balls of one color obeys 
in the limit a normal distribution. Gouet (1993) further shows, under the 
assumptions of Bagchi and Pal, convergence of the discrete urn evolution 
to a stochastic Gaussian process, and he also investigates other cases using 
martingale arguments. Aldous, Flannery and Palacios (1988) observe that 
such results can be supplemented by almost sure convergence properties — 
their treatment extends the relation between branching processes and urn 
models to be found in the book by Athreya and Ney [(1972), Section V.9]. 
Thanks to the works of these and many other authors, the normal evolution 
of the process in the central regime can thus be regarded as well understood. 

In this article, we revisit urn models under the radical angle of analysis. 
[Aldous (1991) otherwise provides an insightful comparison of the scopes of 
the traditional probabilistic approach and the modern methods of analysis 
of algorithms in his introductory section.] Our main results provide a com- 
plete analytic solution describing the composition of the urn at each instant, 
but, although our methods potentially apply to all the 2x2 balanced urn 
schemes, we focus attention in this paper on urns involving subtraction, that 
is, having negative diagonal entries. The matrix can accordingly be taken 
under the form (— a,a + s;b + s,— b), with balance s > and diagonal coef- 
ficients —a, —b < 0. (The urn's initial composition is fixed with to balls in 
total of which ao are of the first type.) Such models with negative diago- 
nal entries are occasionally mentioned by some authors as a harder nut to 
crack, since the direct embedding of urn schemes into branching processes 
explained in the book of Athreya and Ney [(1972), Section V.9] ceases to be 
directly applicable. [This position is perhaps to be taken with caution given 
the discussion in Aldous, Flannery and Palacios (1988) of extensions of the 
classical probabilistic framework.] 

In the first part of the article (Sections 1 and 2), we introduce the partial 
differential equation approach to urn models with two types of balls and 
constant row sum. Our analysis starts with a partial differential equation 
(PDE) that is linear of the first order and that describes exactly snapshots 
of the urn compositions at all times. The solution of this partial differential 
equation, obtained by the standard method of "characteristics," provides 
an indirect expression for a bivariate generating function that encodes the 
possible configurations of the urn at each time n. It is found that this bivari- 
ate generating function is expressed in terms of a fundamental function tp, 
which is defined implicitly by an equation of the form 

(2) i>(I(u)) = Q(u,v). 

There (u, v) lies on a Fermat curve u h + v h = 1 with h = a + b + s a sort of 
"complexity index," the quantity Q being a rational function on the curve, 
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and I(u) an Abelian integral on that same curve — that is, the integral of 
a rational function on the curve. The parameterization (2) suffices in all 
cases to determine the dominant singularities of ip together with the as- 
sociated singular expansions. As a consequence, analytic principles provide 
the (known) Gaussian law for the urn's composition at large instants, to- 
gether with a precise determination of the speed of convergence as well as 
an explicit form of the large deviation function in terms of the Abelian in- 
tegral I(u) [see (3) for a specific instance]. In general, ip is associated with 
algebraic curves of genus strictly higher than 1. (Note: The Fermat curve 
is of high topological genus [Lang (1982)], namely g = (h — l)(h — 2)/2, so 
that one already has g = 10 in the case of the 72,3 model discussed below. 
This makes the occasional existence of elliptic function solutions, which are 
objects of genus 1, quite remarkable.) 

Our investigations were initially motivated by a desire to understand the 
specific urn model 7^3 := _g), which forms the subject of the second 
part of this article. This particular urn process intervenes as a model of 
several schemes for managing an important data structure of computer sci- 
ence known as the search tree [Knuth (1998) and Mahmoud (1992)] and 
it surfaces in the analysis of 2-3 trees and fringe-balanced binary search 
trees [Aldous (1991), Aldous, Flannery and Palacios (1988), Bagchi and Pal 
(1985), Eisenbarth et al. (1982), Panholzer and Prodinger (1998) and Yao 
(1978)]. What is striking about this urn is that the model can be completely 
resolved in terms of elliptic functions of the Weierstrafi type. For instance, 
our general results express that the probability of large deviations at time n 
is exponentially small in n with a rate that is a simple transform of the 
integral 



A parallel elliptic connection had been uncovered earlier by Panholzer and Prodinger 
(1998) using rather different methods. Their penetrating analysis depends on 
the specific relationship that the 7^3 model entertains with a special type of 
"fringe-balanced" search trees — a root decomposition of the tree then leads 
to a perturbed nonlinear ordinary differential equation (of the rough form 
Y'" = Y' 2 + • • • ) akin to the one satisfied by the Weierstrafi p-function. In 
this particular case, our elliptic connection for the 7^3 urn model could al- 
ternatively be deduced by reverse-engineering of the Panholzer-Prodinger 
treatment, combined with an easy reduction of a special urn model studied 
by Mahmoud (1998). We do not proceed along those lines since Panholzer 
and Prodinger's nonlinear differential approach is problem-specific, and, for 
example, it would not yield the other elliptic cases listed in (4). 

The general character of our analytic results, Theorems 1 and 2, actually 
permits us to single out all the cases where elliptic solutions prevail, namely, 
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all urns of balance 1 or 2 and a sporadic urn of balance 3, corresponding to 
the six matrices 



(4) 
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We have chosen to illustrate the specificity of the elliptic models in a con- 
crete way by developing properties of the 72,3 model (Section 3) and then 
proving our classification theorem regarding the elliptic cases (Theorem 4 
and Section 4). 

1. Analytic solution of the general case. We now take up the general 
case of a balanced urn model with two types of balls and negative diagonal 
entries. The matrix is of the form 

(5) M=(- + a s °_ + /), a,„> , 

with s > the balance. Start with ao balls of the first type ("black") and 
bo balls of the second type ("white"), so that to = a o + bo is the initial size; 
the size of the urn at time n is then exactly the deterministic quantity 
tn = to + ns. In order for the urn not to be blocked by an infeasible request, 
the usual "tenability" conditions [Bagchi and Pal (1985) and Gouet (1993)] 
for urns with subtraction are assumed: 

,„s J (To) : a divides ao and b divides bo; 

\ (Ti) : a divides b + s and b divides a + s. 

We shall see soon that all such models are "solvable by quadrature" in the 
sense of Taylor [(1996), page 86]. In other words, only elementary algebraic 
functions, composition and inversion, as well as integration are involved in 
the solution, as is expressed by the general statement of Theorem 1. There 
results a complete characterization of dominant singularities, as summarized 
by Theorem 2. Probabilistic consequences are subsequently explored in Sec- 
tion 2. 



1.1. Algebraic approach. Based on formal operator calculus, there is an 
elegant symbolic approach to the derivation of PDE's for urn models, which 
establishes a transparent connection between the combinatorial structure of 
a model and the PDE that expresses it. 

The combinatorial model considers all balls involved in the game to be 
distinguished by distinct integer stamps: balls present at time are stamped, 
say, 1, . . . , ao for type B and ao + 1, • • • , to for type W. New balls are stamped 
with "new" numbers: the balls that are taken away from the urn are (conven- 
tionally) the ball selected as well as others taken according to a deterministic 
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policy, for example, by starting from smallest numbers. For instance, the urn 
( - 3 2) initialized with two balls of type B stamped with 1 and 2 may give 
rise to an evolution history starting as 



Time 


1 


2 


3 


choose 2 


choose 3 


choose 6 


choose 1 


Urn 1b, 2b, 


1b,3w,4w, 


, * , 

1b, 5b, 6b, 7b, 


, * , 

1b,5b,7b,8w,9w, "' 



with subscripts indicating colors/types of the corresponding balls. 

In what follows, we consistently use [z n ]f to denote the coefficient of z n 
in the formal power series or analytic function /. 

One first needs to relate combinatorics and probability. We let X n be 
the random number of balls of type B at time n, and denote by p n (u) its 
probability generating function (PGF). Let h n (u) be the counting generating 
function of the evolution histories of length n, where u marks the number 
of balls of type B: the coefficient [u k ]h n (u) is the number of histories com- 
prising n transformations of the urn and resulting in k balls of type B. We 
have po(u) = u a ° as well as ho(u) = u a ° , and in general 

t (t + s) ■ ■ ■ (t + (n - l)s) 
since the total number of possible histories of length n is 

(8) t (t + s) ■ ■ ■ (t + (n - l)s) = n\s n ( n + ^ ~ 1 ) , 

as results from multiplication of n elementary choices. (Naturally, the bal- 
ance condition is crucial to this connection.) Introduce finally the exponen- 
tial generating function of the h n (u), so that 

z n 

(9) H(z,u):=Y,h n (u)-r 

n>0 H - 

is a bivariate generating function (BGF). As u — > 1, the bivariate generating 
function H (z, u) degenerates into a simple algebraic function H (z, 1) = (1 — 
sz)~ to / s , since it then only counts histories in accordance with (8). Thus, 
H(z,u) is a priori a "deformation" of a simple algebraic function. 

For u a variable, we let d u = Jjj be the corresponding partial differential 
operator. It is notationally convenient to make use of the modified operator 

df 

9 u = ud u so that O u f = u—. 

ou 

Differential operators are well known to correspond combinatorially to a 
"pointing" operation. For instance, one has 



d u u a = au a - 1 , 6 u u a = au a 
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so that d u may be interpreted as "select a u-element in all possible ways and 
remove it" while 6 U means "select a u-element in all possible ways and keep 
it." There are many instances in the combinatorics literature of such a us- 
age of differential operators; see, for example, Bergeron, Labelle and Leroux 
[(1998), Section 2.1], Flajolet and Sedgewick (2003) and Goulden and Jackson 
[(1983), page 45]. 

Consider now an urn model defined by a matrix M of the form (5), and 
represent momentarily a particular urn configuration with A white balls and 
/i black balls by the monomial rriA,^ = u^v^. The partial differential operator 
(associated to M), 

(10) r = u - a v s+a e u + u s+b v- b e v , 

is such that the application of T to rriA,^ describes all the possible successors 
of the urn represented by vax^ = u^v^ when one step of ball replacement is 
performed. 

Start with an urn of initial type (ao, bo) represented by u a °v b ° . Let h n (u, v) 
be the polynomial describing all possible evolutions of the urn in n steps. 
[In particular, h n (u,l) = h n (u).] Then, one has 

h n (u,v) = T n o( y u ao v bo ). 
We opt for exponential generating functions and define 

^ -,71 

H(z,u,v) = V h n (u,v)—. 

n>0 

One has symbolically 

H(z,u,v)=e zT o(u a °v bo ), 
where the exponential of operators is defined in the usual way: 

n>0 

Then, the definition of the exponential immediately implies the differential 
relation 

d z (e zT og) = Te zT og. 
In other words, H satisfies the PDE 

(11) d z H = ToH. 

The last equation is almost the PDE we are looking for but not quite 
(it has a supplementary variable, v). Given the balance condition, the urn 
population increases by exactly s at each step. Accordingly, H involves three 
variables, u,v and z, but their exponents in H are bound by a homogeneity 
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condition, each monomial generated being of the form u^v^z 71 with A + fi = 
sn + to- In other words, each monomial m composing H satisfies 

(12) (9 u + v -s9 z )m = t o m, 

and the relation extends by linearity to H itself. 

In summary, a system of two equations now determines H (with 9 U = ud u ): 

(13) d z H = To H, 
{ ' (e u + d v -s9 z )H = t H. 

One can then eliminate the explicit differential dependency on v (the oper- 
ator d v ) and get from (10) and (13) 

d z H = u~ a v 1+a 8 u H + u l+b v l ~ b {se z H - 6 U H - t H). 

At this stage it becomes possible to set v = 1, that is, completely eliminate 
the redundant variable v itself. In this way one obtains the fundamental 
PDE 

(14) (1 - szu b+s )^- + (u b+s+l - u 1 "*)^ - t u b+s H = 0, 

where H = H(z, u) = H(z; u, 1) . 
The main result is then: 

Theorem 1. Consider the urn specified by matrix a -b)> with initial 
conditions (ao,&o) andto := ao + bo, assumingitto be tenable. The probability 
generating function at time n of the urn's composition is 

F K ' s n T(n + t /s) 1 1 v ' h 
where the bivariate generating function H(z, u) is given by 
H(z, u) = 5{u) t0 ^{z5{u) s + /(«)), 

with 

5{u) := (1 - u h f' h , I{u) := £ dt, h:=a + b + s, 

and the function tp is defined implicitly by 



^(J(u)) 



u 



a-o 



5(u 



to 



Proof. We make use of the classical method of characteristics exposed 
in most textbooks, for example, Zwillinger [(1989), Section 94]. Following 
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this method, one first associates to the linear first-order partial differential 
equation (14) the ordinary differential system 

dz du dw 

(15) 



1 - SZU b+s u s+b+1 - u l ~ a t U b+s W ' 

where w "represents" H , and look for its first integrals. 

The equation binding w and u allows for separation of variables, 

dw u* 1 " 1 

— =*o-ft 7-du, 

w vr — 1 

so that a first integral of (15) is 
(16) w8(u)- to =Ci. 

The equation binding z and u is similar but inhomogeneous: 

dz u h - 1 u a ~ l 



sz—r — - + 



du u h — 1 u h — 1 

The homogeneous equation is solved by separation of variables as z = £ • (1 ■ 
u h )~ s / h . By the variation-of-constant technique, one finds 

z = £(«)(! - u h )-'/ h , = - f ^ ITLm/, 



(1 _ t h}(a+b)/h 



so that a first integral of (15) is 

(17) z5(u) s +I{u) = C 2 . 



According to the method of characteristics, the general solution to the 
fundamental PDE (14) is obtained by coupling the two first integrals (16) 
and (17), namely 

&(H(z, u)5{u)- to , z5(u) s + I(u)) = 0, 

for an arbitrary bivariate function Solving symbolically for H puts the 
solution in the form 

(18) H(z,u) = 5{u) to ip{z5(u) s + /(«)), 

for an arbitrary univariate function tp. The initial condition H(0,u) = u a ° 
finally identifies tp as defined implicitly through inversion of I(u), namely, 
ip(I{u)) =u ao /5(u) t °. 

We observe next that ip{z) is analytic at 0. Indeed the tenability conditions 
of (6) imply that a must divide ao and a must divide b + s, hence a divides 
h = a + b + s. In particular, the general form of the parameterization of 
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ip near is ip(u a ) xu a °, that is, ip(z) x u ao//a , which is compatible with 
analyticity. In fact, the expansions involved are of the form 

\ j>0 ) Y?>0 / 

for some real coefficients Aj,/ij and u ranging in a small enough complex 
neighborhood of 0. Examination of the exponents involved in the inversion 
shows that ip(z) can be expanded as a power series in z, and analyticity of if) 
at results. 

□ 

Sensitivity to initial conditions. When the initial state of the urn is 
changed, the functions involved still live in the same general class. Indeed, 
the ip function corresponding to an initial urn of composition (ao,&o) f a<> 
torizes, in accordance with Theorem 1, as 

(19) ^(z) = Mz) ao/a ^ii(z) bo/b , 
where ipi,ipn are determined implicitly by 

(20) *(,<„)) _(^)', *,(/<-)) -(jij)', 

corresponding to an urn initialized with to = clq = a and to = ^o = b, respec- 
tively. The analytic treatment given below extends to both functions tpj,ipjj, 
and it is seen that the main determinant of the category of special func- 
tions encountered is the index h of the Fermat curve and the integral I(u). 
Equations (19) and (20) thus give us flexibility for the choice of the initial 
conditions, as is done repeatedly below. 

In the case where a and b are each at least — 1 , balls have a "descendance" 
and the evolution of descendants are combinatorially independent. Accord- 
ingly, the factorization (20) can be viewed as expressing the fact that the 
histories of all the initial balls can be freely shuffled. (It is known that shuf- 
fle products correspond to products of exponential generating functions.) A 
parallel decomposition underlies the probabilistic reduction of this class of 
urn models to multitype branching processes [Athreya and Ney (1972)], at 
least in the case where no diagonal entry is below —1, so that the disappear- 
ances of balls are not coupled. 

1.2. Complex-analytic structures. For notational simplicity, we shall adopt 
in this section the initial conditions ao = to = a, that is, the urn is initial- 
ized with exactly a balls of the first type (B): by (19), (20) and the ensuing 
remarks, no essential loss of generality is implied by such a choice. 
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1(0) P =I(l) 

Fig. 1. The elementary kite is the image of a small sector So- 

We make use of the quantity h = a + b + s. The function v = 5(u) corre- 
sponds to the complex Fermat curve, 

u h + v h = l, 

which has topological genus g = (h — l)(/i — 2)/2. Following a classical ter- 
minology, the integral I(u) = J u a ~ 1 v~ a ~ b is a particular Abelian integral 
over this curve. The diagram that summarizes the parameterization of ifj is 
then 

u 



z — > WV Z ) 

The major characteristics of an urn model turn out to be determined by the 
nature of the map u ^ I(u) in the complex plane, with J(u) playing only a 
secondary role. 

As observed in the proof of Theorem 1, the function ip is analytic at 
and it satisfies ip(z) x z a °/ a there. Also, the nature of the parameterization 
near 0, where I(u) x u a implies that I(u) effects an a-fold covering of a 
neighborhood of the origin and that ip(z) is of the form 

(21) i/>(z) =z ao/a 4>{z h ' a ), 

for some ip analytic at the origin. In other words, in order to define ip 
parametrically by means of u, it suffices to let u range in a sector H of 
angle 2n/(h/a) at the origin, and from now on, we shall do so. (As already 
noted, the tenability conditions precisely imply that a divides h.) 

Consider first the complex plane with h rays emanating from and having 
directions given by all the hlh. roots of unity. The sector Sj is defined as 

S j :=(z,z = Re ie ,0<R<oo,^<9<^±^ 

h h 
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We claim (and prove below) that the image of Sq by I(u) is the interior of 
a quadrilateral (Figure 1), with vertices at the points 

0, 1(1), /(+oo), I(e 2 ^ h ), 

and call this quadrilateral /C the elementary kite. [Note: The incomplete 
Beta integrals that make up tp are related to hypergeometric functions as 
well as to the Schwarz-Christoffel integrals of conformal mapping theory. 
For the latter aspects, see, e.g., the book of Nehari (1975), his Exercise 4, 
page 196, and his Chapter V.] One has 



o (l-t h )( a + b )/ h h \h'hj hT((a + 8)/h)' 

where use has been made of the usual Eulerian Beta integral [Whittaker and Watson 
(1927)] 

(22) B(a,(3):= [\ a -\l - tf~ l dt - T( - a ™ 



T(a + f3) ' 

We henceforth denote the quantity 1(1) by p. 

The local mapping properties corresponding to the four vertices of the 
elementary kite are determined by the local behavior of I(u): (i) at 0, I(u) 
multiplies angles by a, so that the angle of the kite at is (ii) at 1, 
I(u) multiplies angles by f , so that the angle of the kite at vertex 1(1) is 
?Y [and similarly for vertex I(e 2ln ^ h )]; (iii) at infinity, I(u) multiplies angles 
by b, so that the angle at I(+oo) is i n order to see that I(u) maps the 
boundary of So to that of /C, observe first that I(u) maps [0,1] onto the 
segment [0, 1(1)} by monotonicity of the integrand. Then, as u continues to 
increase along (1, +oo) passing above 1, the function 5(u) becomes a complex 
number of fixed argument ^^vr. In other words, I maps the ray [l,+oo] to 
the segment [7(1), 7(+oo)] (with +oo here understood as lying inside So). 
A similar discussion gives the mapping properties associated with the other 
two sides of the kite K,. 

Next, we turn to sectors S\, Let £ := e 2i7r//fe . The image of sector Sj is 

simply obtained as the image of So by I(uC J ), which, by a linear change of 
variables [since I(uC J ) = C,~' ia I( u )\ is the image of the elementary kite under 
a rotation of angle see Figure 2 for a particular instance. Because 

of (21) and the accompanying remarks, it is sufficient to consider < j < -. 

Definition 1 . The fundamental polygon of an urn model is the (closure 
of) the union of h/a regularly rotated versions of the elementary kite about 
the origin. 



We state: 
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Fig. 2. The elementary kite (in black) and the fundamental polygon associated with the 
urn (—1,4; 4, -1). 

Theorem 2. Consider a balanced 2x2 urn with subtraction as in The- 
orem 1 and let it be initialized with a$ = a, bo = 0. The corresponding func- 
tion ip is analytic for z in the fundamental polygon of Definition 1. Further- 
more, it is analytic in \z\ < p, where 




(l_^)(a+6)/ft h {h'hj hT({a + s)/h)' 



On \z\ = p, the function ijj is singular at p and at the points puj J where 
u> = exp(2z7r^) is an (h/a)th root of unity, regular at the other points. Its 
singular expansion as z — > p is of the form 

(23) ip(z) = s- t0 ' s {p - z)~ t0 / s A((p - z) h/s ), 

with A analytic at 0, .A(O) = 1, -A'(O) ^ 0. (Principal determinations as 
z — > p~ are assumed.) This expansion extends to a sector of opening larger 
than tv at p. 

At the points z = poJ 3 , the singular expansion is determined from the ex- 
pansion at z = p by the fact that i^(z)z^ aa / a is invariant under the mapping 

Z I — > UJZ. 

The expansion (23) gives ip(z) as the product of a main singular part 
of the form (p — z)~ to / s multiplied by a Puiseux series, that is, a series in 
fractional powers of (p — z). We shall occasionally refer to the quantity h/s 
as the Puiseux exponent of ijj. It plays a special role in the discussion of 
elliptic urns in Sections 3 and 4, in which case it reduces to an integer 
value. 

Proof of Theorem 2. First, the fact that I(u) assumes each value 
in /C once and only once when u £ Sq is a consequence of basic properties of 
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conformal mapping theory, which we recall. Let (3 be an arbitrary number 
interior to K,. The number v({3) of times that I(u) assumes the value (3 6 JC 
for u interior to So is by the residue theorem 



IdSo I{u) - P 

where dX represents the boundary of a region X oriented positively. Then, 
the change of variables I(u) = x gives 

1 f dx 



2in JdK x - (3 

where the reduction to the value 1 is due to the fact that (3 is by assumption 
interior to K,. This implies that the functional inverse u = I ly ~ l \z) is well 
defined (and analytic) for z interior to /C, and so is Tp(z) since tp(z) = J(u) 
while J(u) depends analytically on u. These properties extend in turn to 
the fundamental polygon by rotations of the base sector. 

We next examine the behavior of ip near p = 1(1), corresponding to u in 
the vicinity of 1 (say, u — > 1~ , to fix ideas). The expansion can be constructed 
by means of a local uniformizing parameter, here, 1 — u = r . Write 

5(y) = A(y)(l-y) l /\ 

so that A(y) is analytic at y = l. By the change of variables u \— > 1 — r h , one 
finds 

1(1) - I(u) = h J\l - y h r l A(l - y h )- a - b y s - 1 dy, 
1 rli/h.s (, , (h(b-a + 2)-a-b)s h 



(24) (i_ T M«o 
J(u) 



s v ' V 2h(s + l) 



A(l-T h Yo 



where now r — > corresponds to u — > 1 (the series expansions proceed by 
powers of T h ). Thus the parameterization is of the form 

p - z = -(h l l h r) s U(T h ), Mz) = (h l ' h Ty tQ V(T h ), 
s 

where U, V are analytic at and U(0) = V(0) = 1. By analytic inversion, this 
shows that there exists a full expansion of the type (23), with A analytic 
at 0. In other words, the point p is a singularity of ip that is a branch point 
with dominant singular exponent equal to —to/h. 

By rotational symmetry, an expansion of a nature similar to (23) also 
holds at the conjugate points puji where u = e 2i?r / /l is an /ith root of unity. 
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Since ip(z) has nonnegative coefficients, it satisfies Pringsheim's theorem and 
is thus analytic in \z\ < p. By the triangle inequality, we have \I(ue ie )\ < I(u) 
for u £ (0, 1) and 9 € (— 7r,7r). Since the nonzero terms composing the Taylor 
expansion of / at the origin are of the form u a+: ' h , the inequality \I(ue td )\ < 
I(u) is strict as soon as 6 is not a multiple of ^ and I(u) is invertible. From 
there, it results that ip is analytic on \z\ = p except for the regularly spaced 
singularities quoted in the statement. 

This provides the analytic continuation of ip in the fundamental polygon 
as well as in the disc of radius p. If the fundamental polygon is such that 
s/h> i , analytic continuation of if) outside its disc of convergence is granted 
and the proof is completed — this is the situation exemplified by Figure 2. 
Otherwise, the convergent character of A provides the analytic continuation 
of tp in sectors rooted at singularities and extending beyond the disc z < p. 
(This last situation is encountered in the T23 model detailed below; see 
Figures 3 and 4.) □ 

For instance, the urn ("^ _f) gives rise to the fundamental polygon dis- 
played in Figure 2. One has s = 3, h = 5 and 5(u) = (1 — u 5 ) 1 / 5 , so that 
the fundamental polygon is a star with five branches. At the origin, we find 

I(u) = u + ygii 6 H and ip(z) = z + j^z G H . There is an algebraic branch 

point at p where tp(p — x) x (p — x) -1 / 3 and at the conjugate points puji 
where = 1. The nature of the branch point of tp at p is 

Mz) = (3Z)~ 1 /3 (1 _ |j(3Z) 5 / 3 - ^(3Z) 10 / 3 + •••), Z:=(p- z), 

the Puiseux exponent associated with A being the fractional number h/s = 
5/3. 

2. Probabilistic consequences. Singularity analysis [Flajolet and Odlyzko 
(1990) and Odlyzko (1995)] makes it possible to extract very precise informa- 
tion on coefficients of a generating function once the function is recognized 
to have isolated singularities on the boundary of its disc of convergence. 
Under such conditions, assuming first unicity of the dominant singularity <r, 
an asymptotic estimate for a function F of the form 

F(z) ~ c(l-z/a)- a 

z — >c 

valid in a complex region beyond a (a sector centered at a of opening angle 
larger than ir and including the disc of convergence) entails a matching 
estimate for the function's coefficients: 

\z n }F{z) ~ ca- n — -. 
1 J v J n^oo r(a) 
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Full asymptotic expansions can be transferred from functions to coefficients 
in a similar way [Flajolet and Odlyzko (1990)]. Also, in the presence of sev- 
eral dominant singularities, contributions to coefficients are to be composed 
additively. (This technology based on Hankel contours is of the complex 
Tauberian type.) It applies to the function ip(z) itself, and since ip(z) = 
H(z,0), it provides immediately sharp estimates of the probabilities that all 
balls are of the same color at epoch n, which corresponds to extreme large 
deviations. 

Corollary 1 (Extreme large deviations). For any balanced 2x2 urn 
with subtraction (i.e., negative diagonal entries), the probability that balls at 
time n are all of the same color and of the second type (W) is 

Proof. The singularity at z = p of H(z, 0) = ip(z) contributes to [z n ]?/>(2;) 
a term 

jto/s— i / / i 



We have the periodicity expressed by (21). Thus, for n in a suitable con- 
gruence class, there are h/a similarly behaving singularities to be combined. 
The total number of histories of length n is, from (8), asymptotic to 

n t /s-l 

n\s n — — r^. 
r(t /s) 

The result follows after normalization by the latter quantity. □ 

Next, we summarize the basic technology used to derive a Gaussian limit 
by the following statement, a simplified form of what is often referred to 
as the "quasi-powers theorem" and originates in works of Bender (1973) 
and Hwang (1998). Throughout this article, we use E and V to denote the 
expectation and variance operators. 

Lemma 1 (Quasi-powers theorem). Let q n (u) = ¥,(u Yn ) be a family of 
probability-generating functions relative to discrete random variables Y n . As- 
sume that there exist two functions A(u),B(u) analytic in a neighborhood V 
of u= 1, such that, in this neighborhood the quasi-power approximation 

(25) q n (u) = A(u)B(u) n (l + £ n (u)) asn^oo 

holds, where \e n (u)\ = 0(n~ 1 / 2 ) uniformly with respect to u, that is, sup ugV \e n 
0(n~ 1 / 2 ). Assume also the variability condition 

VY 

(26) <xV where a 2 := lim — — = B"(l) + B'(l) - B'il) 2 . 

n— >oo n 
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[Equivalence between the two forms of a 2 is granted under condition (25).] 
Then, the random variables Y n converge in law to a Gaussian limit, with 
speed of convergence 0(n~ 1 / 2 ): for any x, one has 

Y n -E(Y n ) _\_ 1 f* ^ 2/2 ^, , 1 



Proof [Sketch; see Bender (1973) and Hwang (1998) for details]. The 
characteristic function q n (e lt ) of the Y n is by assumption closely approxi- 
mated by an nth power. The variable Y n is next centered around its mean 
and scaled by its standard deviation in the usual way. A calculation similar 
to the usual case of independent random variables [e.g., Billingsley (1986), 
page 367] then shows the standardized version of q n {e lt ) to converge to 
e~* I 2 , which is the characteristic function of a Gaussian law. The speed of 
convergence estimate finally results from the Berry-Esseen inequality found 
in Lukacs (1970). □ 

In order to apply the quasi-powers theorem, we choose a small complex 
neighborhood of u = 1 and keep u in this neighborhood. The BGF H{z,u) 
rewrites as 

(28) H(z,u) = 5(u) to ip{p-5(u) s (K(u) - z)), 
where 

1 r 1 t" 1 ^ 1 

( 29 ) K ( U ) : =TT^ F7 \ , dt, 
v 1 y ' S(u) s Ju 5(u) a+b 

and K{u) has a removable singularity at 1 with K{\) = 1/s. Treating u 
as a parameter, we find that, as a function of z, the quantity H(z,u) has a 
singularity at z = K(u) that gets smoothly displaced when u varies. Because 
of the nature of the singularity of tjj at p, the singular exponent remains equal 
to the constant — to/s. Thus, for some function L(u) that is analytic at u = 1, 
one has 

[z n ]H(z, u) = L{u)K{u)- n n t0 / s ~ l (l + O (-^ 

the error term being uniform by virtue of uniformity of the singularity anal- 
ysis process [Flajolet and Odlyzko (1990)]. This has the shape of a bona fide 
quasi-powers approximation for the probability generating function, 



[z n ]H(z,l) L(1)\K(1)J V \n h / s 
The quasi-powers theorem then applies and gives: 
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Corollary 2 (Gaussian law and speed). For any balanced 2x2 urn 
with subtraction, the random variable X n representing the number of balls 
of the first color (B) at time n is asymptotically Gaussian with speed of 
convergence to the limit 0(re -1 / 2 ), as expressed by (27). 

The fact that the limit distribution is Gaussian was first observed by 
Bagchi and Pal (1985). These authors applied the moment method and de- 
termined the main asymptotic orders of moments of the centered variable 
X n — E(X n ). Their method does not, however, appear to give access to the 
speed of convergence as expressed above. This speed is on the other hand 
neatly implied by the functional limit theorem of Gouet (1993). Here, we 
emphasize that the speed of convergence comes out almost immediately from 
the analytic approach. 

In general, the moments are computable systematically from the exact 
expressions of Theorem 1 by successive differentiation with respect to u upon 
setting u = 1 and making use of the singularities of ip and its derivatives as 
expressed by Theorem 2. All moments happen to be expressible in closed 
form. 

Corollary 3 (Moments). For any balanced 2x2 urn with subtraction 
and any r > 0, the rth factorial moment of the distribution of X n is of 
hyper geometric type: it is a finite linear combination of terms of the form 

/n+to / s+t—kh/ s— 1\ 

} , 0<k,£<r. 



^n+to/s— 

The existence of such finite binomial forms for moments of all orders 
does not seem to have been previously noticed. Explicit forms are given by 
Kotz, Mahmoud and Robert (2000), but only for the first moment and at the 
cost of some labor, in the case of the urn model (4,0; 3, 1). Bagchi and Pal 
(1985) obtained such expressions for a wide class of urns, but in the case of 
the first two moments only. 



Xr (z) := -(d r u H(z,u)) u=1 = [(u - iy\H{z,u) 



Proof of Corollary 3. The quantity 
1 

is a generating function of the rth factorial moment of X n in the sense that 

E( * J- ra^)' 

where the notation X- is the usual notation for falling factorials [Graham, Knuth and Patashnik 
(1989)], namely, 

(30) a L = a(a-l)---(a-r + l). 
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In order to gain access to such moments, we make use of the singular 
expansion of ip near p. Given the variant form of H{z,u) from (28), the 
singular expansion (23) of Theorem 2 provides the alternative representation 

(31) H(z, u) = s- to l s {K{u) - zy to/s A((l - u h )(K(u) - z) h/s ), 

which is our starting point. This expansion is analytically valid when (say) 
\z\ < p/2 provided u stays in a small enough neighborhood of 1. Write A{w) = 
J2k>o a k wk ■ One then has 

(32) H(z, «) = - u h ) k {K{u) - z ) kh /s-to/s_ 

k>0 

Clearly, for the rth moment, it suffices to consider the sum in (32) with the 
index k restricted to values in the interval [0, r], so that 

(33) Xr(z) = 1 £ a k (ff u ({l - u h )\K(u) - z)-^ s+kh ' s )) u=1 . 

k=0 

The function K{u) is analytic at u = l. Accordingly, the quantity (K(u) — 
z)^ 1 and its derivatives at u = 1 are of the form 

(Km- r 1 - _L_ ° 2K 'M ^ 3 k'(i) 2 s 2 k"{i) 

{ {) ' ~1-8Z' (l-Sz) 2 ' (I-SZ) 3 (I-SZ) 2 ' 

and so on, with similar formulas holding for fractional powers. Thus Xr(z) 
is invariably an algebraic function of a very special form, namely a finite 
linear combination of terms of the type 

(1 - sz yto/s+kh/s-e^ o<k,£<r. 
The statement then follows by coefficient extraction. □ 

As a consequence, one gets mechanically, 

mtv \ s + k / v \ sh 2 (s + a)(s + b) 

s + h (s + h) z (s + 2h) 

which is consistent with the estimates of Bagchi and Pal [(1985), pages 395- 
397]. 

Finally, we turn to large deviations, for which the book of den Hollander 
(2000) can serve as a smooth introduction. It is known from the works of 
Hwang (1996) that a quasi-power approximation (in the sense of Lemma 1) 
for a family of PGFs leads to very precise "moderate deviation" estimates 
valid in some range not too far from the center of the distribution. We recycle 
here the technology of Hwang (1996), though the range is a little different. 
The large deviation rate is fully characterized by the following statement: 
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Corollary 4 (Large deviations). Consider any balanced 2x2 urn with 

s+h> 



subtraction. Let £ be any number of the open interval (0, s|4t). One has 



(34) lim - log P(X n < £ ■ n) = -K(g) , 

n^oo n 

where the rate function 1Z is determined from K(u) defined in (29) by 

(35) K(£)= max \og(sX^K(X)). 

AG(0,1) 

Equivalently, 

(36) ^(0 = log( S A^(A )) 

where Ao € (0,1) satisfies - '—p- °^ + £ = 0. 

[Put another way, TZ((,) is the Legendre transform of log(si^(e*)).] 

Proof. Notice that E(X n ) ~ f+fisn, so that (34) quantifies the left part 

of the distribution as approximately given by e~ n ^^ . The basic ingredient 
is Cramer's technique of "shifting the mean" conjugated with upper bounds 
of the saddle point (equivalently, Chernoff) type as well as lower bounds 
based on the quasi-powers theorem in a shifted region. 
First, one has 



1 



since multiplication by (1 — u) -1 sums coefficients of generating functions. 
Next, for any f(u) analytic at having nonnegative Taylor coefficients, the 
easy inequality [u k ]f{u) < /(A)A~ fc holds provided the positive quantity A 
is taken inside the disc of convergence of f(u). There results from these two 
observations the majorization 

Pn(A) 



(37) ¥(X n < £n) < 



1 - A)AL«™J ' 



valid for any A £ (0, 1). 

In order to derive an upper bound on large deviations, it suffices to 
choose (as usual) the best possible value of A in (37). Now, for fixed posi- 
tive A G (0, 1), the function H(z, A) has a dominant singularity of the alge- 
braic type at z = K(X), see (23). A simple calculation based on the fact that 
the dominant singularities of tp are at poj 3 and that I{u) increases from to 
1(1) = p for u £ (0, 1) shows further that H(z, A) has for A G (0, 1) a unique 
singularity at if(A) on \z\ = K(X). Therefore one has by straight singularity 
analysis, 

(38) Pn (X) ~ C x ■ s- n K(X)~ n , 

n — >oo 
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for some constant C (depending smoothly on A). When £ lies in any fixed 
compact subinterval of (0, sfr|), the upper bound (37) can then be rewritten 

as 

HX n < in) < Cs- n K{\)- n \-^ n 

for some constant C. This is a form amenable to optimization. Let Ao be such 
that K(X)~ l X~^ attains its minimum over (0, 1) at Ao- General convexity 
properties of probability generating functions imply that Ao exists and is 
unique. 

The value of Ao is obtained by cancelling the derivative of K (A)~ 1 A~£ 
and is thus a root of the second equation in (36). Up to factors that are 
sub exponential in n, the upper bound in (38) is of the form e~ n7 ^), with 
7£(£) as given by (36) and (35). We have thus established "one half" of (34), 
namely, 

-logP(X n <£n)<-ft(£)+ (l), 
n 

with 7£(£) determined by (36). 

There finally remains to argue that the upper bound is tight, that is, 
derive a lower bound on the probability values. This results from Cramer's 
technique of shifting the mean. The shifted law r n ^ = [u k ]r n (u) defined by 
the probability generating function 

r n (u) := 

satisfies a standard quasi-powers approximation and is itself amenable to 
Lemma 1. Assume first that the variability condition (26) holds for the 
shifted law given by r n (u) . In that case the sum of probabilities J2^ n -^/n<k<^n r n,k 
of the shifted law tends to a nonzero constant as it is approximated by a 
Gaussian integral. By construction, the r n ^ are the p n ^ = [u k ]p n (u) multi- 
plied by a quantity Aq which varies between e -0(^) Aq > and 0(l)Af . Thus, 
the corresponding sum J2^ n -^/n<k<^nPn,k is, up to subexponential factors 
(themselves of the form e~°( v/ ™)), of the type e~ n ^^\ This implies a lower 
bound, hence the "other half" of the equality in (34). Finally, if the variabil- 
ity condition at Ao is not satisfied (this can only happen at isolated points), 
then an even stronger type of concentration holds for the shifted distribution 
r n ^] in that case, the variance of the shifted distribution is o(n), which, by 
Chebyshev's inequality, entails the stated lower bound on the sum of the 
r n fc, hence the lower bound on partial sums of the p n k- d 



The dual regime of large deviations on the right tail of the distribution is 
determined upon exchanging the roles of quantities a and b. 
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3. The 12 5 3 "tree" model. The model is determined by its matrix and 
the initial conditions 

T 2 ,3=(^ ^3 J , a = 2, b = 0. 

It has motivated much of the study of subtractive urns over the past two 
decades, given its relevance to several data structures of computer science 
[Aldous, Flannery and Palacios (1988), Bagchi and Pal (1985), Eisenbarth et al. 
(1982), Panholzer and Prodinger (1998) and Yao (1978)]. In this section, we 
arrive at the elliptic connection expressed by Theorem 3 and closely related 
to earlier works of Panholzer and Prodinger (1998). Our reason for treating 
this example in detail is twofold: first, it serves as a concrete illustration of 
the general treatment of Sections 1 and 2; second, it paves the way to our 
eventual characterization of the elliptic urn models in Section 4. 

3.1. Basic analytic structure. Taking ao = 2 and bo = corresponds to 
to = 2. Theorem 1 provides an expression for H(z,u): 

H(z, u) = 5{u) 2 ^{z5{u) + /(«)). 

Here h = 6, so that 

S(u) := (1 - u 6 ) 1/6 , I(u) := T dt = [ U t — n r dt, 



k 5{tf Jo (1-t 6 ) 5 / 6 
and the function tjj is defined implicitly by 



ip(I(u)) = J{u) where J{u) :- 



u 2 u 2 



5(u) 2 (l-u 6 ) 1 ^' 

The results of Section 1 apply directly to this case. The elementary kite, 
which is the image of the sector So of opening ^, is a quadrilateral with 
vertices at (0, p, I(ooe 2t7T / 12 ), puj), where cj := e 2l7r / 3 . We find after a simple 
computation, upon following the proper branch of 5, 

I(+ooe 2 ^ 2 ) = p-e^ f' w2dw 

Jo (1 



i)5/6 



e^ 1 - B ( 1 -, 1 -)=p(l-e^^ 



P ~ 6 V6'27 r V " 2 

Thus, as u varies from to +oo, passing through 1 (and above it), I(u) 
describes first the segment from [0,p], then the segment [p, p(l + u)/2]. The 
kite in this case happens to be a triangle and we shall refer to it as the 
"elementary triangle" (Figure 3). 

There is also a "double parameterization" [due to evenness of both I{u) 
and J(u)], so that we may freely identify points u and —u. To this effect, 
we define 

(39) H := > 0) V ((3(z) = 0) A > 0))}. 




1 1(1) = p 

Fig. 3. The "elementary kite," here a triangle, To (right) is the image of the basic sector 
So (left) via the mapping u \— > I(u). 



Then, as z ranges over Ti, the fundamental polygon is obtained by gluing 
three rotated images of the elementary triangle. It is thus an equilateral 
triangle with center at the origin (Figure 4) — we call it the "fundamental 
triangle." 

Next, the local analysis of ip at its dominant singularity z = p results from 
the general treatment offered in Section 1. We find 

(40) i(j(z) = Z~ 2 -\Z 4 + ^ 7 Z 10 + ---, Z:=p-z. 

What is noteworthy here is the presence of a pole, rather than an algebraic 
singularity that prevails in the general case covered by Theorem 2. Similarly, 
the points puj and puj 2 are double poles, so that the function 

(41) ~ ( TTz? + I^W + (fJ- ) 

is analytic in a disc \z\ < R for some R> p. 

The fact that the dominant singularities of ip are poles naturally led us 
to look for the next layer of singularities, as this would provide very precise 
information on the exponential smallness of error terms. In so doing, much 




Fig. 4. The "fundamental polygon," here a triangle, T (right) is the image of the slit 
upper half-plane (Ti.) (left) via the mapping u>—*I(u). 
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to our surprise, we uncovered a lattice structure commonly associated with 
elliptic functions (Figure 5). 

3.2. The elliptic structure. An elliptic function is a function that is mero- 
morphic in the whole complex plane and is doubly periodic. Amongst the 
many different ways to develop the corresponding theory, perhaps the sim- 
plest is the one originally proposed by Weierstrafi, where elliptic functions 
are defined as sums of rational functions taken over lattices. [Accessible 
introductions appear in the books by Whittaker and Watson (1927) and 
Chandrasekharan (1985).] 

Definition 2. A lattice A with generators £,77 G C is defined as the set 
of complex numbers 

A(£,r?) = {niC + n 2 r?|ni,n 2 G Z}. 
The Weierstrafi p-function relative to A is classically defined as 

(42) p(*;A) = ^+ E ■ ((7^-^ 

(The Weierstrass p-function is by construction doubly periodic.) 

We shall make use here of the "hexagonal" lattice A defined as the lattice 
generated by e 47r / 6 , e~ tn ^ 6 , see Figure 5, 

(43) A hcx := {me i7r/6 + n 2 e-^ /6 |m^2 g Z}, 

and its associated Weierstrafi zeta function, p(z; Ah ex )- We state: 
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Theorem 3 (Elliptic connection). The ^-function of the 7^3 model ini- 
tialized with two balls of the first type (do = to = 2) is exactly 

1 (z~P\ ■ , 1 r(i/3)r(i/6) 
(44) ^wl^bl) " tt " := 6 r ( i/2) ■ 

where p{z) := p(z;Ah ex ) is the Weierstrafl function of the hexagonal lattice. 
In particular, the bivariate generating function of the model is expressible in 
terms of elliptic functions. 

Note. The function if) can be alternatively written as ip(z) = p(z — 
p | 0, —4) where p is specified by the lattice invariants gi = and g% = —4. 

Proof of Theorem 3. Consider the whole complex plane tiled by 
nonover lapping copies of the hexagon of center p, radius pV3, having vertices 
at the points p + p^A^. 

We claim that any complex point z is reachable as a value I(j(u)), where 
the notation 7(7(7/)) indicates that the integral defining / is to be taken 
along a path 7(11) that starts at and ends at u. Similarly, J(^y(u)) will 
represent the determination of J(u) along path 7(7/) that is obtained by 
continuity from the principal determination at 0. Otherwise said, we are 
walking on the Riemann surface of the Fermat curve 5(u). 

The algorithm is as follows. Assume for simplicity that z is the center of 
one of the equilateral triangles in which the hexagonal tiling decomposes. 
The straight line Lq from to z can be first slightly deformed into a curve L\ 
that avoids all the vertices of the tiling. This L\ can then be transformed into 
a polygonal line L2 that connects centers of successive equilateral triangles. 
Finally, each segment of L2 can be changed into a pair of segments going 
through one of the vertices of the lattice and forming an angle a positive 
multiple of ir/3. The resulting polygonal line, L3, will be called the standard 
z-path. See Figure 6 for a graphic rendering. 

The contour 7, called the standard u-path, is then obtained from the stan- 
dard z-path L3 by first applying a contraction by a factor p, then executing 
the following routine: 

• turn by an angle of 69 whenever L3 turns at an angle of 9 (where 9 is a 
multiple of tt/3) around a vertex of the lattice, 

• turn by an angle 9/2 (where 9 is a multiple of 27r/3) whenever L3 turns 
by 9 around the center of one of the equilateral triangles. 

The construction is then easily modified to accommodate points that are 
not centers of triangles of the tiling. 

For any z in C that is not a vertex of the tiling, the algorithm described 
above determines constructively a path j(u). By design, along such a path, 
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Fig. 6. A standard path in the z -plane from to P = z and the contour 7 above the 
u-plane that realizes it via u t— ► z = 1(7(11)). 

one has I(7(ii)) = z. Indeed, the standard u-path is precisely such that it 
"undoes" the effect of I(u) on angles at points either vertices of the tiling 
or centers of the triangles; at the same time, the variation of I(u) along a 
segment from a point no above to a point u\ above some C J is precisely 
of modulus p and thus gives rise to a segment with the "right" length. See 
once more Figure 6. In this way, we find that I{^f{u)) reaches any point z of 
the complex plane that is not a vertex of the tiling, and at the final point, 
J(l( u )) 1S locally analytic, so that tp is itself analytic at z. Thus, ip(z) can 
be continued to the complex plane punctured at vertices of the tiling. 

When z = w is one vertex of the lattice, then it is approached from a 
certain direction by a path "f(u), where u is near £°, or £ 2 . Along the 
path a certain determination 5°(u) of 5 is in force, where all determinations 
are of the form ( r 5(u) with < r < 6. Then, the very same determination 6° 
must be adopted in J{^(u)) that tends to infinity as I{^{u)) approaches w. A 
local analysis entirely analogous to the one conducted for the three dominant 
poles shows that ip has a double pole at w, and that its principal part there 
consistently exhibits the same dominant coefficient and residue. 

The analytic continuation of ip{z) along such paths 7 therefore has the 
same dominant parts and residues at double poles as the right-hand side of 
(44), namely the function 



entire function reduces to results from Liouville's theorem, as we finally 
argue. 

Draw discs of some sufficiently small but fixed radius around the six roots 
of unity in the u-plane and consider these as excluded regions in the con- 




Consequently, the difference 




That this 
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struction of li-paths. Then the image z = I(j(u)), as j(u) varies, avoids the 
plan stripped of small ovals around the corresponding lattice points of the 
z-plane. But, at the same time J(^(u)) remains bounded. Therefore, on the 
complex plane with "holes," ip(z) is uniformly bounded by a constant. From 
the fact that ip is doubly periodic, there results that it is also bounded over 
the plane with holes, hence 



for some c\ > 0. In particular, the bound holds on an infinity of near-circular 
contours centered at the origin and having arbitrarily large diameter. Then, 
by virtue of a known variant of Liouville's theorem (an entire function 
bounded in modulus along large contours is a constant), one must have 
identically 



for some complex constant d\. This constant is^actually equal to as is 
seen from comparing the expansions of i/j(z) and ip{z) at 0. The proof of the 
theorem is completed. □ 

3.3. Probabilistic consequences of the analytic model for T2.3. We are 
now in a position to exploit the analytic solutions expressed by Theorem 3. 
The general theory of Sections 1 and 2 applies, giving the large deviation 
rate function and the limit Gaussian law. In addition, curious exact repre- 
sentations as sums over lattice points result for the probability generating 
functions describing the urn composition (Section 3.3.1). Surprisingly per- 
haps, a very precise form of all moments can be obtained in terms of a family 
of polynomials of "binomial type" [Rota (1975)]; see Section 3.3.2. 

3.3.1. Exact representations and Gaussian laws. The lattice structure 
that underlies the Weierstrafi function is directly reflected at the level of 
coefficients. The resulting form below is naturally very strong, as it is an 
exact description of the probability generating function at time n. 

Corollary 5 (Elliptic structure of T23). For the 72,3 model, the proba- 
bility generating function p n (u) = K(u Xn ) admits an exact formula valid for 
all n> 1, 



\ip(z) - tp(z)\ <c x 



ip(z) -rp(z) =di 




where 
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Proof. From Theorem 3, we need to extract [z n ]5 2 ijj{5z + I), where 
ip admits a decomposition as a sum of rational fractions over elements of 
the lattice A. Then, after a simple calculation, one gets 



n + *-.f?i*\ {w - z) z j 



u>eA 

where A* is a translated and scaled version of A 



A- = ^A + A». 



The result follows. □ 



Corollary 2 regarding general Gaussian limits applies here. Thus, for the 
7^3 model, the random variable X n representing the number of balls of the 
first type at time n is asymptotically Gaussian with speed of convergence 
to the limit 0(n -1 / 2 ), in the sense of (27). The random variable X n superfi- 
cially resembles a sum of independent random variables since its probability 
generating function is essentially an nth power of the fixed function K{u)~ l . 
It is, however, of interest to observe that the function K{u)~ l , though ana- 
lytic at and satisfying K (1) = 1, is not a probability generating function, 
as its Taylor coefficients of index 6, 12, 18, . . . turn out to be negative: 

K{u)~ l = 0.713 + 0.254u 2 + 0.090u 4 - 0.086u 6 + 0.022u 8 + • • • . 



3.3.2. The shape of moments. An interesting consequence of the ellip- 
tic connection concerns moments of the distribution of the urn's composi- 
tion, X n . Bagchi and Pal, Mahmoud, and Panholzer and Prodinger have 
determined the exact form of the first two moments, while Bagchi and Pal 
have obtained further asymptotic information on the moments of higher 
order. This already involved a certain amount of calculational effort with 
recurrences. In fact, globally, the moments have an amazingly simple form 
deriving from the elliptic connection. 



Corollary 6 (Moments of 72, 3)- For the 7^3 model, exact polynomial 
forms of moments of any order are available: the factorial moments satisfy 

E((X n )^)=P r (n + 2), n>&r-l, 

where the P r are polynomials generated by 
00 , r 

(46) e vL{h) =Y, -\Pr(v) and L(h) = - log If (1 + h). 

r=0 r ' 
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Using a symbolic manipulation system, the polynomials are easily com- 
puted from the expansion of K at u = 1. To wit: 

K{l + h) = l-^h + f l h^^-^ + .... 

One then finds mechanically 

PM = —i P 2 H = — (521/ + 17), 

Paiy) = ^y( 1976z/2 + 1938u ~ 11063 )- 

In particular, the mean and variance of X n are 

E(X n ) = |(n + 2), V (X„) = |§(n + 2) 2 . 

Proof of Corollary 6. Take the fundamental PDE, isolate G' u (z,u) 
and repeatedly differentiate with respect to u, then set u = 1. This provides 
a triangular system from which one can "pump" in succession the generating 

functions of moments of order 1,2,3, One then verifies by induction that 

the ordinary generating function of the moments of order r is of the form 

where P r ,Q r are polynomials and 

deg(P r (z)) <r, deg(Q r (z)) < 6r - 2. 

This argument grants us nonconstructively the existence of a polynomial 
representation for each moment as soon as n is large enough. 

There remains to identify the particular class of polynomials involved. 
Start from the fact that 

p n (u) = K(u)~ n ~ 2 + exponentially small terms in n. 

Since the factorial moment of order r satisfies 

E(X£) = (d r uPn (u)) u=1 = [(« - 1)»), 

it can be obtained, up to exponentially small error terms, by expanding 
K{u)~ n ~ 2 around u=l. Retaining only the polynomial part (in n), 

[(« - l) r )K{u)- n - 2 = [{U - l)r] e -(n+2)logX(u) = ^ e ~(n+2)logK(l+h) ^ 

we get what the statement asserts. □ 
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3.3.3. Large deviations. Corollary 4 applies to the effect that the large 
deviation rate function is a transform of K{u). An immediate consequence 
of the analysis of the polar singularities of ip is a very precise quantification 
of extreme large deviations: 

Corollary 7 (Extreme large deviations of ^2,3)- The probability that, 
at time 3n + 1 in the 7^3 model, all balls are of the second color (W) is (any 
A<8): 

3p- 3n - 3 (l + 0(A- n )). 

Proof. The function ip{z) is exactly the BGF of the urn at u = 0. Thus 
(n + l)" 1 [z n ]-i/>(z) is the probability for the urn not to contain any ball of 
the first type. The property then results immediately from the fact that the 
next layer of poles is on \z\ = 2p. □ 

4. The elliptic cases. In this section, we list all the cases of urns with 
subtraction which, like the 72,3 model, lead to elliptic functions; such urns 
are, as in the previous section, attached to exact lattice representations and 
simplified moment forms. We shall say that an urn is elliptic if, for some 
choice of initial conditions (ao, bo), the fundamental function ip is a fractional 
power of an elliptic function, 

(47) ^(z)=n(s;A)*/9, 

where II is meromorphic and admits a period lattice A. The power will 
depend on the initial conditions, but by taking the initial configuration de- 
termined by ao and &o sufficiently large, one can always render the exponent 
integral; see the remarks on "sensitivity to initial conditions" as well as 
(19) and (20). Obviously, for urns with subtraction, one need only consider 
models that are arithmetically irreducible: by this is meant that the matrix 
M = (" $) defining the urn has coprime entries: gcd(a, (3, 7, 5) = 1. 
The key characters in this section are the following six urns: 

- 3? S 3? si: 3? il 

The urn A is of course the 7-2 3 model. As is easily checked, any urn of balance 
s = 1 is necessarily of type A, B or C and any arithmetically irreducible urn 
of balance 2 can only be of type D or E. Thus, the first five cases exhaust all 
possible types of urns with subtraction having balance s = 1,2. The urn F 
is one of the four possible irreducible urns having balance s = 3. We state: 
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Fig. 7. The six elliptic cases in order A,B,C, D, E, F : the diagrams formed by the fun- 
damental polygon together with its rotated images. (The elementary kite is darkened.) 

Theorem 4. All balanced 2x2 urns with balance s = 1 (cases A, B and 
C) are elliptic. All urns with balance s = 2 (cases D,E) are also elliptic. 
There exists one 11 sporadic" urn (case F) of balance s = 3 that is elliptic. 
These six matrices represent the only arithmetically irreducible models of 
urns with subtraction that are elliptic. 

Proof. For 2-3 trees, the matrix A corresponds to a = 2,6 = 3,s = 
l,h = 6. As seen in Section 3.2, it is associated to a regular tiling of the 
plane by equilateral triangles. Its ip function has, by virtue of Theorems 2 
and 3, a double pole at p and all other points of the lattice and is a Weierstrafi 
p- function. 

Next, we turn to the other urns of balance 1. The corresponding mapping 
properties are represented in Figure 7. The model B has an elementary 
kite which is a right triangle with vertices 0,p,ip, so that the fundamental 
polygon is the square with vertices p,ip, —p, —ip. This fundamental square 
tiles the plane. The model C leads to a lattice much similar to the 2-3 
tree case. Theorem 2 shows that in all the three cases A,B and C, the 
ip function has a pole at p since the principal exponent at the singularity 
to/s as well as the Puiseux exponent h/s of (23) are integers (here s = 1). 
Double periodicity then results from arguments similar to those developed 
in the proof of Theorem 3. 

An analogous discussion applies to the two urns of balance 2, namely D 
and E. In the case of D, the function I(u) is even exactly a lemniscatic 
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integral: 

v ' Jo 

Finally, among the four urns of balance 3, namely 



-1 4 \ f-1 4 \ f-1 4 \ f-2 5 

4 -lj> \ 5 -2J' \ 7 -4 J ' ^8 -5 

one found to be elliptic has s = 3, h = 6, so that /i/s = 2 is integral: this 
is i 7 = ( _ g _ 2 ). By the usual reasoning, the function ^ leads to an elliptic 
function when one starts with oq = 3 and bo = 0. 

Note a necessary condition for an urn to be elliptic: in addition to the 
tenability conditions, s should divide h, in accordance with Theorem 2 and 
(23); that is, the Puiseux exponent h/s should be integral. (Otherwise, all 
powers of tp inherently have branch points and hence cannot be meromor- 
phic.) The other three cases of balance s = 3, including the "pentagonal" urn 
of Figure 2, correspond to a fractional value of the Puiseux exponent h/s in 
(23) and therefore cannot be reduced to elliptic functions. 

There finally remains to prove that all elliptic urns have indeed been 
found. The condition that the Puiseux exponent h/s is integral, arithmetic 
irreducibility, and the tenability conditions taken together imply the exis- 
tence of triples (x,y,z) of integers, representing up to possible permutation 
the values (a,b,s), such that 

(49) gcd(x, y, z) = 1, x \ y + z, y \ z + x, z \ x + y. 

Simple arithmetic shows that the only values for which the system admits 
a solution are permutations of the basic types 

(50) (1,1,1), (1,1,2), (1,2,3), 

and, in particular, one must have s < 3. The arithmetic argument goes as 
follows. Take x <y < z. One has z < 2y since z < x + y by the third divisibil- 
ity condition in (49); then note the stronger property that x,y,z are pairwise 
coprime (proof: a contrario). Then set z + x = qy, where q < 3; the first di- 
visibility condition then implies x\{q + l)y, and since x and y are coprime, 
x|4 so that x S {1, 2, 4}. This in turn implies y < z < y + 4 while z must di- 
vide y + x. Combining this with the second and third divisibility conditions 
of (49), we see that there are only finitely many possibilities which are then 
easily tested. Finally, completeness of the list (50) implies that elliptic cases 
have indeed all been found. □ 



It is pleasant to note that the elliptic urns correspond to the crystallo- 
graphic groups of the Euclidean plane, that is, groups of isometries acting 
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discontinuously (in fact, the ones which admit a compact fundamental do- 
main). As is well known, these groups themselves describe the possible reg- 
ular tessellations of the plane by polygonal tiles and are in finite number; 
see, for example, Berger (1977) and Yoshida (1997). 

5. Discussion. A probabilist may legitimately expect more than standard- 
issue central limit theorems and Cramer approximations to come out of an 
analytic treatment of urn models. We would like to offer a few comments. 

1. The combinatorial derivation to the fundamental PDE of (14) applies 
(with only notational adjustments) to any urn model with more than two 
colors, provided it remains balanced. The resulting PDE is invariably first 
order and linear. This study has exhibited a class, the subtractive balanced 
2x2 models, for which the associated ordinary differential system provided 
by the method of characteristics and the corresponding generating functions 
prove to be analytically tractable. This suggests to look for other cases where 
analysis can be made to work. 

2. On the algebraic-analytic front, Theorems 1 and 2, though they appear 
to admit no universal r x r extension when r > 2, can at least be adapted 
and generalized to several models [work in progress with V. Puyhaubert 
(2003-2004)]: 

(i) 2x2 balanced urns with positive entries of the type of "Bernard 
Friedman's urn" [Friedman (1949)]. Interestingly enough, an adapted form 
of the algebraic solution expressed by Theorem 1 holds (only integration 
constants need to be changed). This may be seen as an elicitation of some of 
Friedman's remarks concerning his differential recurrences [Friedman (1949), 
pages 61 and 62]. Developments parallel to Theorem 2 appear to be possible, 
they then open access to non-Gaussian laws which have not been previously 
made explicit. 

(ii) 2x2 balanced urns corresponding to nonnegative triangular matri- 
ces. Gouet (1993) shows strong functional limits to exist, but some of the 
involved characteristics have remained inaccessible due to nonconstructive 
aspects of martingale theory. In that case, we can supplement Gouet's find- 
ing and derive explicitly stable laws and a local limit theorem. 

(iii) Balanced urns with three colors and a nonnegative triangular matrix, 
as well as some special cases of triangular r x r urns for r > 3. 

3. The determination of the large deviation rate is obtained by standard 
arguments of Cramer type, once the relevant analytic forms have been es- 
tablished, but, to the best of our knowledge, it is new. We observe that 
the asymptotic approximations we have obtained, when combined with the 
saddle point method, can provide strong forms of large deviation estimates, 
complete with subexponential factors and multipliers. As we have seen in 
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the case of moment estimates and Gaussian laws, a precise determination 
of the speed of convergence to the asymptotic regime is a boon provided by 
the analytic machinery. We also observe that the general algebraic solutions 
supplied by Theorem 1, when coupled with (19) and (20) (sensitivity to ini- 
tial conditions), make it possible, in principle, to analyze the evolution of 
the urn starting with a large number of balls of each color. This could be 
put to use in order to analyze finely the distribution of sample paths, in the 
central or large deviation regime. 

4. Besides integer partitions and theta functions, elliptic models occasion- 
ally pop up in combinatorics, some models related to permutations having 
been found by Dumont, Flajolet, Frangon and Viennot about 1980. The 
present work contributes a new kind, urn histories, themselves equivalent 
to certain weighted lattice paths. On another register, given that algebraic 
functions of higher genus can be uniformized by theta- Fuchsian functions, 
we may even fantasize about the possibility of "nonprobabilistic" repre- 
sentations for urns of genus higher than 1 in terms of tessellations of the 
hyperbolic plane. 
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